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Abstract. Thermalization play a central role in out-of-equilibrium physics of 
ultracold atoms or electronic transport phenomena. On the other hand, entanglement 
concepts have proven to be extremely useful to investigate quantum phases of matter. 
Here, it is argued that bipartite entanglement measures provide key information on 
out-of-equilibrium states and might therefore offer stringent thermalization criteria. 
This is illustrated by considering a global quench in an (extended) XXZ spin- 1/2 
chain across its (zero-temperature) quantum critical point. A non-local bipartition 
of the chain preserving translation symmetry is proposed. The time-evolution after 
the quench of the reduced density matrix of the half-system is computed and its 
associated (time-dependent) entanglement spectrum is analyzed. Generically, the 
corresponding entanglement entropy quickly reaches a "plateau" after a short transient 
regime. However, in the case of the integrable XXZ chain, the low-energy entanglement 
spectrum still reveals strong time-fluctuations. In addition, its infinite-time average 
shows strong deviations from the spectrum of a Boltzmann thermal density matrix. In 
contrast, when the intcgrability of the model is broken (by small next-nearest neighbor 
couplings) , the entanglement spectra of the time-average and thermal density matrices 
become remarkably similar. 
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1. Introduction 

Rapid progress in the field of ultracold atoms [T| offer brend new perspectives to 
realize controlled experimental setup to investigate out-of-equilibrium physics. Real- 
time observation of quantum dynamics of isolated systems have become possible [2]. 
In addition, ultracold atoms loaded on optical lattices P H] or laser-cooled Coulomb 
crystals of charged ions jl] offer very clean experimental implementation of simple lattice 
many-body Hamiltonians and provide simulators for Condensed Matter. The ability to 
dynamically change parameters [5] in these Hamiltonians on short time scales could be 
exploited to realize quantum information processing pQ or cooling [B] devices. 

In electronic condensed matter systems, relaxation towards steady states play a 
central role in many transport phenomena, like e.g. in electric transport resulting from 
the application of a sudden voltage bias at the edges of a quantum dot [7j or of a 
Hubbard chain j8]. Spin chains also offers simple generic systems to investigate out-of- 
equilibrium physics as e.g. heat transport [9J. However, conceptually, thermalization JTU] 
of non-equilibrium isolated quantum many-body systems after e.g. a sudden change of 
Hamiltonian parameters (quantum quench) is still poorly understood, despite recent 
work on correlated bosons [TT1 IT2] in one-dimension (ID). Generally, whether some 
local observables approach steady values and whether their time average equal the 
corresponding thermal average are often used as criteria of thermalization. However, 
the fact that thermalization occur for certain local observables (according to the above 
criteria) does not at all guarantee that other observables will also meet the criteria. 

Independently, quantum information concepts have been applied with great success 
to several domains of Condensed Matter [T3], giving new type of physical insights 
on exotic quantum phases. Quantum entanglement of a A/B bipartition of a many- 
body (isolated) quantum system can be characterized by the groundstate (GS) reduced 
density matrix (RDM) pa obtained from tracing out the B part. The corresponding 
entanglement Von Neumann (VN) entropy SVn = ~^{pa In Pa} offers an extraordinary 
tool [2], e.g. to identify underlying conformal field theory (CFT) structure in one- 
dimensional systems. Another central quantity is the entanglement spectrum (ES) 
defined by the (positive) eigenvalues of a (dimensionless) pseudo-Hamiltonian "H defined 
by the implicit relation pa = exp (— 7-L). Remarkably, ES faithfully reflect CFT 
structures [15], topological symmetries (16] or properties of edge states in fractional 
quantum Hall states [T7| or low-dimensional quantum magnets |18j . 

So far, time evolution of entanglement has been investigated only in very simple 
cases e.g. for a small segment in a ID system after global or local quenches f2H[ [TH| IT3]. 
However, the full potential of entanglement measures has not been fully exploited to 
investigate thermalization of many-body systems. Because for a given bipartition (i.e. 
into two halves) of the whole system all the information of the GS is contained in 
its Schmitt decomposition and the ES is a (convenient) way of arranging the Schmitt 
coefficients, the ES contains the whole information of the state. The main goal here is 
therefore to use the bipartite ES to take the place of local observables to investigate 
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thermalization: whether the ES approach the spectrum of a thermal density matrix 
(whose effective temperature is fixed by the excess energy brought in the quench) 
is therefore proposed as a stringent test of thermalization. If the ES satisfies this 
thermalization criterium, other observables would also do. The proposed choice of a 
bi-partition of the total system is crucial since (i) an extensive subsystem is considered 
so that thermalization of local and non-local observables is addressed at once, (ii) finite 
size scaling can be performed and (iii) reduced density matrices benefit from numerous 
conserved quantities (total momenta and total spin) allowing for a direct comparison of 
the ES separately in each symmetry sector and hence an ultimate comparison between 
time-averaged and Boltzman thermal density matrices. Note that non-local real-space 
partitions have also been used e.g. to define non-local order in gapless spin chains [T9] . 

Practically, a genuine correlated anisotropic ID spin chain and the time evolution 
of its reduced density matrix, entanglement entropy and ES after a global quench are 
considered using exact (Lanczos and full) diagonalization techniques. For the integrable 
case we have considered, the reduced density matrix exhibits strong time-fluctuations 
and its infinite-time average significantly deviates from a thermal density matrix, despite 
the fact that the entanglement entropy reaches a well-defined entropy plateau. In 
contrast, thermalization, as defined by the above criteria, seems to be possible when 
integrability is broken by adding (small) extra terms to the Hamiltonian. 

2. Model and setup 

Let us now consider the ID anisotropic spin- 1/2 Heisenberg (so-called XXZ) model 



whose (ID) parameter space can be mapped on a (half) unit circle assuming J xy = cos9 
and J z = sin 9. Alternatively the system can be viewed as a (half-filled) hard-core 
boson chain with hopping t = J xy /2 and nearest neighbor (NN) repulsion V = J z . Its 
remarkable phase diagram obtained by Haldane |2T] and shown in Fig. [l](b) exhibits a 
Quantum Critical Point (QCP) located exactly at the SU(2)-symmetric point 9 = 7r/4 
separating a gapped Charge Density Wave (CDW) insulating phase (Ising phase in the 
spin language) and a gapless metallic Luttinger Liquid (XY phase). Phase Separation 
(PS) occurring for an attractive interaction \V\ > It (9 < — vr/4) will not be of interest 
here. Ultimately, we shall consider adding next NN hopping t' and repulsion V (see 
Fig. [IJa)) in order to break integrability. 

The key to construct extensive quantities is to realize a non-local partition of the 
chain into "even" and "odd" sites. In other words, if the chain is drawn in a zig-zag 
fashion as in Fig. [IJa), the A and B parts form the two edges of the system which become 
explicit for t' ^ and V ^ 0. One consider here finite chains of N (= 16, 20 and 24) 
sites with periodic boundary conditions, as shown in Fig. [2j The groundstate RDM of 
the subsystems Pa = Pb can be computed using translation symmetry of each L = N/2 
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Figure 1. (Color online) (a) A N-site (extended) XXZ spin-1/2 chain (drawn here as 
a "zig-zag" on a periodic ribbon) is partitioned into two identical A and B subsystems 
of L = N/2 sites by cutting bonds along the dashed line, (b) Phase diagram of the 
XXZ spin (or hardcore boson) chain mapped onto a circle assuming J xy = cos 9 and 
J z = sin#. (c) The CDW groundstate for 9 — ir/2 i.e. t — (J xy — 0) prepared before 
the quench : all bosons (up spins) are located on A and the B sites are empty (down 
spins) . 

site subsystem [18J (each symmetry class is labelled by a momentum K = 2nn/L) and 
the conservation of the Z-component of the total spin (i.e. the number of bosons), 
S\ + S§ = 0. Interestingly, liquid and insulating bosonic phases are characterized by 
qualitatively different entanglement properties. First, the entanglement entropy per 
unit length, shown in Fig. [2ta) for a N = 24 site chain, shows a (kink-like) maximum at 
the SU(2)-symmetric QCP and vanishes (in the thermodynamic limit) for the classical 
CDW (Ising) configuration obtained when 9 — > ir/2. Indeed, as shown in Fig. []Jc), for 
9 = n / 2 the (two-fold degenerate) GS is a simple product of a completely filled (A or 
B) chain by a completely empty (B or A) one. Note that the symmetrized GS (i.e. the 
equal weight superposition of the two classical GS) still retains a 1/L finite size entropy 
as shown in Fig. [2]^a). Secondly, each quantum phase is uniquely characterized by its 
ES defined as the spectrum of "H = — lnp^: Fig. |2](b-c) (Fig. [2^d-e) ) for parameters in 
the CDW phase (LL phase) shows very distinctive features and a clear gapped (linear 
gapless) spectrum. In addition, at the QCP (Fig. [2^d)) one observes a SU(2) multiplet 
structure. 
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Figure 2. (Color online) GS properties of the XXZ chain; (a) VN entanglement 
entropy vs 9 computed on a N=24 site ring. The entropy is normalized by the maximum 
value Smax = Lin 2 (L — N/2), (b-e) Typical entanglement excitation spectra (for the 
4 values of shown in (b)) as a function of momentum K along the ribbon. The 
eigenvalues £ Q of T~L = — In pa are labelled according to the Z-component of the total 
spin (i.e. the number of bosons Na — L/2 — Sa) of the A subsystem (legend of symbols 
on graph). Note the CDW (LL ), doubly-degenerate (unique) GS of H (of energy £o) 
carries Na = or Na = L (Na = L/2) hardcore bosons. 
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Figure 3. (Color online) (a-c) Squarcd-fidelity (shaded) and VN entanglement entropy 
of various non-equilibrium states obtained after different sudden quenches of the 
N = 24 sites chain Hamiltonian, 9i nit — > Of as shown on plots (t' = V = 0). The 
continuous lines (dashed red line) correspond to a symmetric (non-symmetric) initial 
state (see text). The dotted (green) line is the GS entropy for 9 = Of. 



3. Time evolution and bipartite entanglement entropy 

Let us now consider a sudden quantum quench of the system at time r = (quasi- 
adiabatic quenches will be treated later on). For simplicity, the initial state |0(O)) is 
chosen either (i) as one of the two (zero entropy) degenerate GS at 9 = 9^ = n/2 (i.e. 
for vanishing hopping) shown in Fig. flic), where all hardcore bosons are located on a 
single edge, A or B, or (ii) as a symmetric combination of the two. At times r > 0, 
the boson hopping t (spin-flip term in spin language) is switched on, i.e. the value of 9 
discontinuously jumps at r = to its "final" value 9{ (see Fig. [j](b)). 
The time evolution of the system wavefunction, 

|0(r))=exp(-irfi-(0 f ))|0(O)> (2) 

is easily computed by (arbitrary) time steps of 5t (from 0.01 to 0.8) with arbitrary good 
precision [22] (typically better than 10~ 16 ) by Taylor expanding the time evolution 
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operator exp (— iSrH(0{)). Hence, for finite size systems under consideration, an exact 
computation of the time- dependent RDM, 

p A (r)=Tr B |0(r))(0(r)| (3) 

can be done. Results for the squared-fidelity F(t) = | (0(0) |0(r)) | 2 and the entanglement 
entropy, 

S VN (r) = -Tt{ Pa (t) In p A (r)} (4) 

are shown in Fig. [3]for increasingly "strong" quenches corresponding to 9f = 2n/5, n/3 
and 7r/6. After a very short transient regime the squared-fidelity drops sharply and 
the entanglement entropy raises to a more or less well defined plateau, whose average 
value exceeds the value of the GS of the final Hamiltonian (f-GS). In the N — > oo limit, 
one expects a non-perturbative regime [23] where F{t) ~ 2~ N — > 0. However, on finite 
system and for small quench, F(r) can remain large, as seen e.g. in Fig. [3]^a). In that 
regime, integrability of the model can play an important role [21]. Therefore, from 
now on, one will assume a sufficiently large quench, let say Of = 7r/6, to observe time 
evolutions on finite clusters generic of the thermodynamic limit. This corresponds in 
fact to a quench "across" the QCP at 9 = tt/A into the region of stability of the LL 
phase. However, from Fig. [3|c) we can see a well defined entropy plateau, suggesting 
that the system does reach a steady state after the quench. Note also that one finds 
that non-symmetric and symmetrized initial states give very similar results so that, for 
simplicity, one can restrict to the symmetric case hereafter. 

The time evolution of any (static) observable O of the A subsystem is given by 
0(r) = Ti(pa(t)0). Time-average like O = (O), where (G) = lim^oo i G(r)dr, 
can then be rewritten as Tr(p^ e O), involving the time-averaged RDM p A e = (pa)- Next, 
time-averages will be performed in a At = 40 time interval using a mesh of 50 points 
(excluding the initial short transient regime). 

Bipartite entanglement can provide precise characterization of the system after the 
quench, in the thermodynamic limit, e.g. showing whether or not it reaches a (quasi-) 
steady state. One finds that the time fluctuations of SVn(t) shown in Fig. ^ a) vanish in 
the thermodynamic limit as revealed by the finite size scaling of Fig. |4j^c) . Interestingly, 
as seen in Fig. Hfa), we notice that the VN entropy of the average density matrix 
p A ve , = — Tr{p^ e lnp^ e }, differs from (SVn), the time-average of the entanglement 
entropy. A proper finite size scaling of the two quantities shown in Fig. |4](b) proves that 
this fundamental property holds in the thermodynamic limit. 

Incidentally, it is also important to distinguish between "quantum" fluctuations 
(O - O) 2 = Ti{p A vc (0 - O) 2 } and "time" fluctuations ((O - O) 2 ). An extreme example 
is the case of S A where S A (r) = at all times (from the time-conserved A <H- B 
symmetry) while quantum fluctuations are finite as shown in Fig. |4Fc). 
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Figure 4. (Color online) (a) Time-dependent VN entanglement entropy after a sudden 
#init = n/2 — > 8f = tc/6 quench (t' = V' = 0), for different chain lengths N as shown 
on plots. Dashed-dotted (dashed) lines corresponds to (Svn) (^vn)- O 3 ) Finite- 
size scaling of the averages shown in (a), (c) Finite-size scaling of the relative time 
fluctuations of the entanglement entropy (red squares) and quantum fluctuations of 
S\ (normalized by L/2). Time-averages are all performed in a At — 40 time interval 
using a mesh of 50 points. The data of b) and c) are well fitted assuming ~ 2~ aL finite 
size corrections. 



4. Entanglement spectra: time evolution within the "entropy plateau" 

In order to understand the physical origin of the difference between (Svn) and Sy^ it is 
interesting to inspect more closely the behavior of Pa{ t ) within the "entropy plateau" 
regime as a function of time r. For this purpose, it is convenient to use the (bipartite) 
time- dependent ES defined by the (positive) eigenvalues of the (dimensionless) pseudo- 
Hamiltonian "H(r) defined as, 

p A (T)=exp(-H(T)). (5) 

As for the equilibrium GS, the spectrum of H. is computed separately in every sector of 
the momentum K = (using translation invariance of the A and B half-systems) and 
of the z-component of the total spin. Fig. [5] shows "typical" entanglement spectra taken 
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Figure 5. (Color online) "Snapshots" at different times r of the ES of pa- (a) lies 
within the transcient regime and (b-f ) lie within the " entropy plateau" . Symbols are 
similar to Fig. [2jb-e) . 



at different times. At short times, some memory of the initial state is clearly visible as 
in Fig. |5^a). At longer times, let's say r > 2.5, it is remarkable to find very different 
spectra at different times although all such spectra lead roughly to a very comparable 
value of the entanglement entropy. In other word, Pa{t) fluctuates strongly in time 
on a constant entropy "hyper-surface" of the space of density matrices with Trp = 1. 
This naturally explains why (SVn) and Sy£ differ substantially. Our finite size scaling 
suggests that this is a fundamental phenomenon and not a finite size effect. 

5. Entanglement spectra: comparison between infinite-time average and 
thermal ensembles 

The RDM pa contains all relevant information about the subsystem, much beyond 
any local observable. As argued in the introduction, its associated ES is therefore 
an ultimate observable to investigate thermalization. Furthermore, the existence of 
conserved quantities such as momentum and particle number makes the comparison 
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Figure 6. (Color online) Entanglement spectra of the time- averaged RDM obtained 
after a 9i n a — tt/2 — > 9f = ir/6 quench in XXZ chains of different lengths, N = 16, 
20 and 24 (i' = V' = 0). Symbols are similar to Fig. [2^b-e) and time-averages are 
performed as for Fig. [4] 



between time-average and thermal density matrices much finer, providing a stringent 
thermalization criterion (based on the close similarity between these two quantities). 

The ES of Pa° computed for XXZ chains (in the entropy plateau) with 16, 20 and 
24 sites, and shown in Fig. [6^a-c), reveal a smooth convergence with system size, an 
accumulation of low pseudo-energy levels at K = and a small pseudo-energy gap, in 
sharp contrast with the ES of the f-GS shown in Fig. [2^e). 

In order to probe thermalization, we now wish to compare p A ve to the thermal 
average of pA- Thermal averages can be defined as: 

^ = E^a|*«><*«l, (6) 

a 

where the prime means the sum is restricted to eigenstates of H(0f) with non- 
zero overlap with |0(O)), accounting for conserved quantities like momentum and z- 
component of total spin, and hidden conservation laws in the integrable case. Both 
canonical (u^ an = ^exp(— (3E a ) where E a are eigenenergies associated to \^ a )) and 
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Figure 7. (Color online) ES of p% c (a,d), p<% n (b,e) and //f cro (c.f) obtained by 
/wZZ diagonalization of ./V = 20 site integrable (t' = V' = 0) and non-integrable chains 
(ftnit = tt/2 -> 6> f = tt/6 quench). Symbols are similar to Fig. [jjfb-e). 1//3 ~ 2.07 (b), 
A£ = 0.3 (c), 1//3 ~ 1.48 (e) and AE = 0.1 (f). 



microcanonical (u;™ lcro = Cst in an energy window [E — AE,E + AE]) thermal 
ensembles can be considered. Note, the effective temperature 1/(3 of the canonical 
ensemble is (implicitly) given by constraining the conserved mean-energy to be E = 
— (V — V')L/2. Note that a full diagonalization of H(9f) is now required since the 
previous formula ^ involves the complete set of eigenstates \^ a ), hence limiting the 
available chain lengths to iV = 20. The infinite-time average p^ e can be computed 
exactly also using ^ with u^ ve = |(\1/ Q |0(O))| 2 which, in contrast to the thermal 
averages, depends now on the initial state. Incidentally, this enables to control the 
very good accuracy of the previous approximate averaging procedure, as shown by a 
direct comparison between Fig. |6](b) and Fig. [7|a). 

The results shown in Figs. [7] reveal striking differences between integrable and 
non-integrable Hamiltonians. For the latter, surprisingly good agreement between the 
low-energy ES of p^ e and p c ^ n is seen, suggesting that the Eigenstate Thermalization 
Hypothesis [lOj may apply up to the level of the (bipartite) entanglement properties of 
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eigenstates. In contrast, for the integrable XXZ chain, noticeable differences between 
the various averages reveal incomplete thermalization which seems to be a fundamental 
property of such systems. 

6. Quasi-adiabatic quenches 

The system after the quench possesses an excess of entanglement entropy (per unit 
length) compared to the f-GS. There is in fact an interesting correlation between this 
excess entropy and the "speed" at which the quench is performed. A smooth quench can 
be realized by considering a time-dependent Hamiltonian H(9(t)) where 9{r) decreases 
continuously from # init = 7r/2 to 9f = tt/6 in a time-interval [0,Tf]. One can choose e.g. 



and T\ = Tf/10. As shown in Fig. [8j under increasing the characteristic time T\ the 
average value of the entropy plateau decreases towards the GS value, as expected 
for increasingly adiabatic processes. Interestingly, the deviation (SVn) — SyN also 
simultaneously decreases (to reach zero for a fully adiabatic process). 

Sudden and continuous quenches also give very different ES of the time-averaged 
reduced density matrix p^ e - Data for continuous and sudden quenches are compared 
in Fig. [9} Remarkably, the spectra in the quasi-adiabatic case resemble very closely the 
equilibrium f-GS spectrum of Fig. 2(e) of the main paper with a clear linear envelope, 
in contast to the spectra obtained for a sudden quench. However, we note that the 
requirement (e.g. on the time scale 7\) to have a quasi-adiabatic process becomes more 
and more stringent for increasing system size because of the vanishing of the finite-size 
gap (in the XY phase). 

7. Summary 

To summarize, the concept of bipartite entanglement is introduced in a genuine out- 
of-equilibrium isolated many-body system. I propose that it provides a simple and 
complete probe of thermalization, beyond the investigation of simple local observables. 
In this framework, absence of thermalization is diagnosed whenever the reduced density 
matrix deviates (once averaged over a sufficiently large time interval) noticeably from 
the thermal (reduced) density matrix taken at an effective temperature imposed only 
by the initial conditions. Generically, we find that the system approaches quickly an 
entanglement entropy plateau. However, integrable and non-integrable chains behave 
quite differently within this plateau. The RDM of the integrable chain strongly 
fluctuates in time and does not thermalize according to the above criterium. This is to 
be contrasted to the case where extra terms are introduced in the chain Hamiltonian 
to break the integrability: in such a case, the ES of the RDM and the spectrum of 




(7) 
(8) 
(9) 




w T = l/(l + exp((T f /2-r)/T 1 )), 
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Figure 8. (Color online) VN entanglement entropy of various non-equilibrium states 
obtained after a 9^ = ir/2 — > Of = ir/6 quench realized on a 24-site XXZ chain 
(t' = V' = 0) over different time scales characterized by T\. The arrows indicate the 
difference between the time- averaged entropy (Svn) (dot-dashed lines) and the VN 
entropy of p^°, S"vn (dashed lines). Time-averages have been performed using 50 bins 
in a At = 40 time-interval within the plateaux. 

the thermal density matrix become very similar, a strong hint that thermalization can 
occur up to the level of an extensive subsystem even though the full system is completely 
isolated. 
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